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Abstract 

We review our previous work on the dynamic structure factor S(k,uj) of liquid Ge (£-Ge) at 
temperature T = 1250 K, and of amorphous Ge (a-Ge) at T = 300 K, using ab initio molecular 
dynamics [Phys. Rev. B67, 104205 (2003)]. The electronic energy is computed using density- 
functional theory, primarily in the generalized gradient approximation, together with a plane wave 
representation of the wave functions and ultra-soft pseudopotentials. We use a 64-atom cell with 
periodic boundary conditions, and calculate averages over runs of up to about 16 ps. The calculated 
liquid S(k, to) agrees qualitatively with that obtained by Hosokawa et al, using inelastic X-ray 
scattering. In a-Ge, we find that the calculated S(k,uj) is in qualitative agreement with that 
obtained experimentally by Maley et al. Our results suggest that the ab initio approach is sufficient 
to allow approximate calculations of S(k, iv) in both liquid and amorphous materials. 
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I. INTRODUCTION 



Ge is a well-known semiconductor in its solid phase, but becomes metallic in its liquid 
phase. Liquid Ge (£-Ge) has, near its melting point, an electrical conductivity characteristic 
of a reasonably good metal (~ 1.6 x 10 _4 n _1 cm _1 yJ), but it retains some residual structural 
features of the solid semiconductor. For example, the static structure factor S(k), has a 
shoulder on the high-fc side of its first (principal) peak, which is believed to be due to a 
residual tetrahedral short-range order. This shoulder is absent in more conventional liquid 
metals such as Na or Al, which have a more close-packed structure in the liquid state and a 
shoulderless first peak in the structure factor. Similarly, the bond-angle distribution function 
just above melting is believed to have peaks at two angles, one near 60° and characteristic 
of close packing, and one near 108°, indicative of tetrahedral short range order. This latter 
peak rapidly disappears with increasing temperature in the liquid state. 

These striking properties of £-Ge have been studied theoretically by several groups. Their 
methods fall into two broad classes: empirical and first-principles. A typical empirical 
calculation is that of Yu et a/j^l, who calculate the structural properties of £-Ge assuming 
that the interatomic potentials in £-Ge are a sum of two-body and three-body potentials 
of the form proposed by Stillinger and Weber Q]. These authors find, in agreement with 
experiment, that there is a high-k shoulder on the first peak of S(k) just above melting, which 
fades away with increasing temperature. However, since in this model all the potential energy 
is described by a sum of two-body and three-body interactions, the interatomic forces are 
probably stronger, and the ionic diffusion coefficient is correspondingly smaller, than their 
actual values. 

In the second approach, the electronic degrees of freedom are taken explicitly into ac- 
count. If the electron- ion interaction is sufficiently weak, it can be treated by linear response 
theoryj^]. In linear response, the total energy in a given ionic configuration is a term which 
is independent of the ionic arrangement, plus a sum of two-body ion- ion effective inter- 
actions. These interactions typically do not give the bond-angle-dependent forces which 
are present in the experiments, unless the calculations are carried to third order in the 
electron-ion pseudopotentialji], or unless electronic fluctuation forces are included jj] Such 
interactions are, however, included in the so-called ah initio approach, in which the forces on 
the ions are calculated from first principles, using the Hellman-Feynman theorem together 
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with density-functional theory^,!?],!^] to treat the energy of the inhomogeneous electron gas. 
This approach not only correctly gives the bond-angle-dependent ion-ion interactions, but 
also, when combined with standard molecular dynamics techniques, provides a good account 
of the electronic properties and such dynamical ionic properties as the ionic self-diffusion 
coefficients. 

This combined approach, usually known as ab initio molecular dynamics, was pioneered 
by Car and Parrinellojj]], and, in somewhat different form, has been applied to a wide range 



of liquid metals and alloys, including £-Ge 
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....... _ , £-Ga r Gei_ r 1 131 , stoichiometric III- 



V materials such as £-GaAs, £-GaP, and £-InP[ 

^-CdTeQ], and ^-ZnTe^sJ, among other materials which are semiconducting in their solid 
phases. It has been employed to calculate a wide range of properties of these materials, in- 
cluding the static structure factor, bond-angle distribution function, single-particle electronic 
density of states, d. c. and a. c. electrical conductivity, and ionic self-diffusion coefficient. 
The calculations generally agree quite well with available experiment. 

A similar ab initio approach has also been applied extensively to a variety of amorphous 
semiconductors, usually obtained by quenching an equilibrated liquid state from the melt. 
For example, Car, Parrinello and their collaborators have used their own ab initio approach 
(based on treating the Fourier components of the electronic wave functions as fictitious 



classical variables) to obtain many structural and electronic properties of amorphous Si 19, 
Isoj ] . A similar approach has been used by Lee and Chang j^jj. Kresse and Hafner 1CJ] 
obtained both S(k) and g(r), as well as many electronic properties, of a-Ge, using an ab 
initio approach similar to that used here, in which the forces are obtained directly from the 
Hellmann-Feynman theorem and no use is made of fictitious dynamical variables for the 
electrons, as in the Car-Parrinello approach. A similar calculation for a-Si has been carried 
out by Cooper et al 22j , also making use of a plane wave basis and treating the electron 
density functional in the generalized gradient approximation (GGA)!^. More recently, a 
number of calculations for a-Si and other amorphous semiconductors have been carried out 
by Sa llk e y et M, and by Drabold and col laboratore Q. These calculations use ai , m Uo 
molecular dynamics and electronic density functional theory, but in a localized basis. A 
recent study, in which S(k) and g(r) were computed for several ab initio structural models 
of a-Si, has been carried out by Alvarez et al[26}. 

Finally, we mention a third approach, intermediate between empirical and ab initio molec- 
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ular dynamics, and generally known as tight-binding molecular dynamics. In this approach, 
the electronic part of the total energy is described using a general tight-binding Hamiltonian 
for the band electrons. The hopping matrix elements depend on separation between the ions, 
and additional terms are included to account for the various Coulomb energies in a consis- 
tent way. The parameters can be fitted to ab initio calculations, and forces on the ions can 
be derived from the separation-dependence of the hopping matrix elements. This approach 
has been used, e. g., to treat l-Sij^. a-Sij3], and liquid compound semiconductors such 
as £-GaAs and £-GaSb[2^]. Results are in quite good agreement with experiment. 

In this paper, we review the application of ab initio molecular dynamics to another 
dynamical property of the ions: the dynamical structure factor, denoted S(k,u). While no 
fundamentally new theory is required to calculate S(k, u), this quantity provides additional 
information about the time- dependent ionic response beyond what can be extracted from 
other quantities [3^. The work we review is the first to calculate S(k,u) using ab initio 
molecular dynamics. Here, we will describe calculations of S(k, u) for £-Ge, where some 
recent expe^Q provided comparison, and also for amorphous Ge (a-Ge). In 
the latter case, using a series of approximations described below, we are able to infer the 
vibrational density of states of as-quenched a-Ge near temperature T = 300K in reasonable 
agreement with experiment. The calculated S(k,u>) for the liquid also agrees quite well with 
experiment, especially considering the computational uncertainties inherent in an ab initio 
simulation with its necessarily small number of atoms and limited time intervals. 

The remainder of this paper is organized as follows. A brief review of the theory and 
calculational method is given in Section II and III respectively. The results are reviewed in 
Section IV, followed by a discussion in Section V. 



II. FINITE-TEMPERATURE DENSITY-FUNCTIONAL THEORY 

The original Hohenberg-Kohn theorem of density-functional theory (DFT) J(| was gen- 
eralized to finite temperatures T e i by Mermin Q]. In this section, we give a brief review of 
the finite-temperature density-functional theory, in a form where it can be used in ab initio 
molecular dynamics calculations. We use atomic units here and throughout this article. 

In the Born-Oppenheimer approximation, the total potential energy U of a system of 
N classical ions, enclosed in a volume V, and interacting with A e ;=NZ valence electrons 
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through an electron- ion pseudopotential v e i(r), is a functional of the valence electron density 
p(r), and also a function of the ion positions (1=1,2,. ..,N). This potential energy can be 
expressed as the sum of a direct Coulombic interaction between the ions and the free energy 
F e i of the electronic system subject to the external potential created by the ions, which we 
denote V ext (r, {Ri}) = Z)i=i v ei(\ r — R-/D- This sum may be written 

U\p(r), {Ri}] = F d [p(r)} + £ W(Ri - RO, (1) 

Kl' 

where iy(R — R') = Z 2 /\H — R'| is the Coulomb interaction between ions of charge Z 
located at R and R'. 

The free energy functional of the electronic system F e i[p(r)] can be partitioned into various 



terms as follows:! 32 
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F el [p(v)} = A s [p(r)] + J p{v)V ext {v)dv + lj P^^l dvdv' + F xc [p(r)]. (2) 

Here the different terms are A s [p(r)], the Helmholz free energy functional for a reference 
system of noninteracting electrons; the interaction energy between the valence electrons 
and the external potential provided by the instantaneous ionic configuration; the classical 
electron-electron potential energy; and the exchange-correlation correction F xc [p(r)} to the 
free energy. A s [p(r)] has contributions from both the kinetic energy T s [p(r)] and the entropy 
tS's [jo(r) ] of the reference system, and may be written 

A s [p(r)]=T s [p^)]-T el S s [p(r)], (3) 

where T e i is the temperature of the electronic system. 

For a given set of ionic positions, the equilibrium electronic number density p(r) is ob- 
tained from the following variational principle: 

J-(U[p(r), {R}] -ft J p(r)dr) = 0. (4) 

Here p is the electron chemical potential, which is chosen to give the desired number of 
valence electrons N d . The Euler equation associated with eq. (jlj) is 

f i = V As (r;\ P }) + V eff (v;[p}). (5) 

Here 



is a functional derivative of and is the "Helmholz free energy potential" (or the kinetic 

potential for T e \ = 0) j^j]. 14// (r; [p]) is an effective one-body potential, given by 

K//(r; [p]) = V^(r) + / ^dr' + ^M. (7) 

It is usual to make several approximations for the functionals A s [p] and F xc [p], and the 
function V ex t(r). First, F xc [p] is only a minor contribution to F e i[p\. As a result, this func- 
tional is often calculated using the local density approximation (LDA) and or the generalized 
gradient approximation (GGA); both are widely used in many applications 32]. 

Furthermore, F xc differs negligibly from its zero-temperature value, the exchange- 
correlation energy E xc , provided T e i < 0.15T/, where Tf is the Fermi temperature [35]. 
Most applications involve temperatures satisfying this inequality. The external field V ext (r) 
contains the electron-ion pseudopotential v e i(r), on which other approximations are made. 
For example, it is common to approximate v e i as local and energy- independent. 

Once an accurately approximated A s [p] is available, one can, in principle, obtain 
Va s (i"; [p]) by functional differentiation, and from this functional, solve eq. (0) to obtain the 
equilibrium electronic density. However, this approach has proved to be very challenging in 
practice. 

Typically, the kinetic energy T s [p] is much larger than the entropy contribution T e iS s [p] 
in the range of temperature interested. Therefore, A s [p] is generally well approximated by 
T s [p\. However, simple local density approximations for T s [p] based on the Thomas- Fermi 
TF) theory are not very accurate, and make several qualitatively incorrect predictions 
3(1 37 1 . More recent developments based on so-called "orbital-free" (OF) approximations 
or the kinetic energy functional T s [p] and the kinetic potential Vr s (r; [p]) 

42 , 3] are still insufficient to make quantitative predictions. On the other hand, although 
these OF methods are currently unable to make accurate predictions, they are computation- 
ally appealing, because the computational time required to use them is independent of the 
number of valence electrons. Therefore, further refinement of such OF methods in DFT may 
lead to valuable new approaches for investigating the properties of large systems, especially 
metallic systems. 

Instead of using approximate forms for A s [p] or Vk s (r; [p]), Kohn and Sham suggested 
a different way to compute the numerical value of A s exactly, by introducing a set of N e i 
orbitals satisfying the coupled KS equations that describe the model system |7[ . Using these 
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orbitals one can determine both A s and the valence electron density p(r). The orbitals ipi 
satisfy the Schrodinger-like equation 

1. 



--V 2 + V eff (r) 



tpi = e^i, (8) 



and are chosen to be orthonormal. The set of equations is made self-consistent by requiring 
that 

p(r) = E/*hMr)-| 2 (9) 

i 

Here fa is the Fermi occupation number, given by 

fi = /(* - A*) = tep(|^) + I)" 1 - (10) 
The kinetic energy T s is given by 

r.[{ifc(r)},{/,}] = £/< / M-^Mdr = £*/*- / p(r)V e// (r)dr, (11) 

and the electronic entropy by 

S.[{/ f }] = -2A; B X^Zn/, + (! - - fi)] (12) 

i 

Eqs. (|7I8I9I1U|) have to be solved self-consistently. The resulting valence electron density can 
be used to computed the electronic free energy F e i[p}; the electronic kinetic energy T s and 
entropy S s are given by eq(JTTJ) and eq (fT2|) 32]. 

The use of finite-temperature DFT is important, in order to eliminate the instability of 
the KS equations at T e i = 0. In practice, T e i is usually regarded as a fictitious tempera- 
ture (chosen for maintaining the stability), and does not necessarily have to equal the ionic 
temperature T. Furthermore, the broadening function in eq. (JIUj) need not be the Fermi func- 
tion, but may be chosen to be some other convenient function. Several possible broadening 
functions have been proposed Q|- 

The forces on the individual ions can now be obtained using the Hellmann-Feynman 
theorem. According to this theorem, the force Fj acting on the I th ion is found b y d iffer- 
entiating the potential energy of the system with respect to the ionic coordinate Rf[45, 46]. 



Specifically, 
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Given the force, one simply applies Newton's second law to get the equations of motion of 
the N classical ions: 

M 7 R 7 = Fj (14) 

In practice, the ions are treated classically, even though the origin of the force is very much 
quantum-mechanical. 

In solving these equations of motion, it is important to maintain the ionic temperature 
T at a fixed value. This is typically accomplished by rescaling the total kinetic energy of 
the ions Ki=Ya=\ qm ^° sa tisfy 

every few time steps. The factor of N-l, not N, appears here so as to conserve the total 
momentum of the system of ions. 

In summary, for each instantaneous configuration of ions, one computes F e i[p(r)} and 
hence U[p(r), {R/}] by solving the finite-temperature Kohn-Sham equations (f7 J) .([5 |) .([0|). and 
(jl(J|) . The Hellmann-Feynman force F/ is then computed from eq. (|13|) . By numerically 
solving Newton's equations of motion for the ions and rescaling the kinetic energy of the 
ions at every time step, one moves the ions to a new configuration. The entire process is 
repeated for a sufficiently large number of time steps to compute the desired property with 
good statistical accuracy. 

III. METHOD 

The method we have used to calculate the dynamic structure factorjscj] is similar to 
that described in several previous papers for other properties of liquid semiconductors jl^. 
Q Q, but uses the Vienna Ab Initio Simulation Package (VASP), whose workings have 
been extensively described in the literature 4?]]." Briefly, the calculation involves two parts. 
First, for a given ionic configuration, the total electronic energy is calculated, using an 
approximate form of the Kohn-Sham free energy density functional, and the force on each 
ion is also calculated, using the Hellmann-Feynman theorem. Second, Newton's equations 
of motion are integrated numerically for the ions, using a suitable time step. The process 
is repeated for as many time steps as are needed to calculate the desired quantity. To hold 
the temperature constant, we use the canonical ensemble with the velocity rescaled at each 
time step. Further details of this approach are given in Ref. 



The VASP code uses the ultrasoft Vanderbilt pseudo-potentials [48], a plane wave basis for 
the wave functions, with the original Monkhorst-Pack (3x3x3) k-space meshes and a 
total of 21,952 plane waves, corresponding to an energy cut-off of 104. 4eV. We have used a 
finite-temperature version of the Kohn-Sham theory 8J in which the electron gas Helmholtz 
free energy is calculated on each time step. This version also broadens the one-electron 
energy levels to make the k-space sums converge more rapidly. The method of Methfessel- 
Paxton of order 1 is used for the broadening function |44j., with the fictitious temperature 
for the electronic subsystem /c s T eZ = 0.2 eV. Most of our calculations were done using the 
generalized gradient approximation (GGA))^ for the exchange-correlation energy (we use 
the particular form of the GGA developed by Perdew and Wangj^), but some are also 
carried out using the local-density approximation (LDA). 

In our iteration of Newton's Laws in liquid Ge (£-Ge), we typically started from the 
diamond structure (at the experimental liquid state density for the temperature of interest), 
then iterated for 901 time steps, each of 10 fs, using the LDA. To obtain S(k) within the 
GGA, we started from the LDA configuration after 601 time steps, then iterated using the 
GGA for an additional 1641 10-fs time steps, or 16.41 ps. We calculate the GGA S(k) by 
averaging over an interval of 13.41 ps within this 16.41 time interval, starting a time t 2 after 
the start of the GGA simulation. We average over all from 1.0 to 3.0 ps. 

For comparison, we have also calculated S(k) within the LDA. This S(k) is obtained by 
averaging over 601 time steps of the 901 time-step LDA simulation. This 601-step interval 
is chosen to start a time t\ after the start of this simulation; the calculated LDA S(k) is also 
averaged over all t\s from 1.0 ps to 3.0 ps. 

To calculate quantities for amorphous Ge (a-Ge), we start with Ge in the diamond struc- 
ture at T = 1600K but at the calculated liquid density for that temperature, as given in 
Ref. Q- 

Next, we quench this sample to 300 K, cooling at a uniform rate over, so as to 
reach 300 K in about 3.25 ps (in 10-fs time steps). Finally, starting from T = 300 K, we 
iterate for a further 897 time steps, each of 10 fs, or 8.97 ps, using the LDA. The LDA S(k) 
is then obtained by averaging over 5.97 of those 8.97 ps, starting a time t\ after the system 
is reached 300K; we also average this S(k) over all t\& from 1.0 to 3.0 ps. To obtain S(k) 
within the GGA, we start the GGA after 5.7 ps of the LDA simulation, and then iterate 
using the GGA for an additional 18.11 ps in 10-fs time steps. The GGA S(k) is obtained 
by averaging over a 15.11 ps time interval of this 18.11 ps run, starting a time t% after the 
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start of the GGA simulation; we also average over all values of ti from 1.0 to 3.0 ps. 

The reader may be concerned that the 3.25 ps quench time is very short, very unrepre- 
sentative of a realistic quench, and very likely to produce numerous defects in the quenched 
structure, which are not typical of the a-Ge studied in experiments. In defense, we note that 
some of these defects may be annealed out in the subsequent relaxation at low temperatures, 
which is carried out before the averages are taken. In addition, the static structure factor 
we obtain agrees well with experiment, and the dynamic structure factor is also consistent 
with experiment, as discussed below. Thus, the quench procedure appears to produce a 
material rather similar to some of those studied experimentally. Finally, we note that exper- 
iments on a-Ge themselves show some variation, depending on the exact method of sample 
preparation. 

We have used the procedure outlined above to calculate various properties of ^-Ge and a- 
Ge. Most of these calculated properties have been described in previous papers, using slightly 
different methods, and therefore will be discussed here only very briefly |50| . However, our 
results for the dynamic structure factor S(k, uj) as a function of wave vector k and frequency 
u) are new, and will be described in detail. We also present our calculated static structure 
factor S(k), which is needed in order to understand the dynamical results. 

S(k) is defined by the relation 

S ( k ) = i<E ex P[* k ■ (*<(*) " r ,(^))])*o - (16) 

y 

where is the position of the i th ion at time t, N is the number of ions in the sample, and the 
triangular brackets denote an average over the sampling time. In all our calculations, we have 
used a cubic cell with iV = 64 and periodic boundary conditions in all three directions. This 
choice of particle number and cell shape is compatible with any possible diamond-structure 
Ge within the computational cell. 

S(k, uj) is defined by the relation (for k ^ 0, w ^ 0) 

1 r°° 

S(k,u) = — y_ oo expM)(p(k,t)p(-k,0))dt, (17) 
where the Fourier component p(k, t) of the number density is defined by 

N 

p(k,t) = ^exp(-zk- ri (t)). (18) 
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In our calculations, the average (...) is computed as 

1 r At i 

(p(k, t)p(-k, 0)> = — J o p(k, U + t)p{-k, t^dU (19) 

over a suitable range of initial times t\. Because of the expected isotropy of the liquid or 
amorphous phase, which should hold in the limit of large N, S(k, u) should be a function 
only of the magnitude k rather than the vector k, as should the structure factor S(k). 

Our calculations are carried out over relatively short times. To reduce statistical errors, 
we therefore first calculate 

S{k, u, h, t 2 ) = -\- t dtp{k, t x + t)p(-k, h) exp(iwt)- (20) 

7TiV JO 

For large enough t 2 , S(k, u>, ti, t 2 ) should become independent of t 2 but will still retain some 
dependence on t±. Therefore, in the liquid, we obtain our calculated dynamic structure 
factor, S ca i c (k, u), by averaging over a suitable range of t\ from to Ati: 

1 /- A 'i 

SaOciKu) = -r— dt 1 S(k,U,t 1 ,t 2 ). (21) 

Ati -Jo 

We choose our initial time in the t\ integral to be 1 ps after the start of the GGA calculation, 
and (in the liquid) At\ = 6 ps. We choose the final MD time t 2 = 16.41ps. For a-Ge, we 
use the same procedure but t 2 = 18. lips in our simulations. For our finite simulational 
sample, the calculated S(k) and 5(k, u) will, in fact, depend on the direction as well as 
the magnitude of k. To suppress this finite-size effect, we average the calculated S*(k) and 
S^k, a;) over all k values of the same length. This averaging considerably reduces statistical 
error in both S(k,u) and S(k). 

Finally, we have also incorporated the experimental resolution functions into our plotted 
values of S(k,cu). Specifically, we generally plot S av (k,u)/S(k), where S av (k,u) is obtained 
from the (already orientationally averaged) S(k, uj) using the formula 

/oo 
R(uu-uj')S(k,uu')duj\ (22) 
-oo 

where the resolution function R{uj) (normalized so that /f^, R{uj)du = 1) is 

R(u) = -^—cxp(-u 2 /lu^). (23) 

In an isotropic liquid, we must have S(k, — u) = S(k,u), since our ions are assumed 
to move classically under the calculated first-principles force. Our orientational averaging 
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procedure guarantees that this will be satisfied identically in our calculations, since for every 
k, we always include — k in the same average. We will nonetheless show results for negative 
uj for clarity, but they do not provide any additional information. 



IV. RESULTS 

A. S(k) for £-Ge and a-Ge 

In Fig. 1, we show the calculated S(k) for £-Ge at T = 1250 K, as obtained using 
the procedure described in Section II. The two calculated curves are obtained using the 
GGA and the LDA for the electronic energy-density functional; they lead to nearly identical 
results. The calculated S(k) shows the well-known characteristics already found in previous 

nn 

simulations j 101. Il2j . Most notably, there is a shoulder on the high-fc side of the principal 
peak, which is believed to arise from residual short-range tetrahedral order persisting into 
the liquid phase just above melting. We also show the experimental results of Waseda et 
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a/|51|; agreement between simulation and experiment is good, and in particular the shoulder 
seen in experiment is also present in both calculated curves (as observed also in previous 
simulations). 

We have also calculated S(k) for a model of amorphous Ge (a-Ge) at 300 A'. We prepared 
our sample of a-Ge as described in the previous section. As for £-Ge, we average the 
calculated S(k) over different k vectors of the same length, as for E-Ge. In Fig. 2, we show 
the calculated S(k) for a-Ge at T = 300, again using both the GGA and the LDA. The 
sample is prepared and the averages obtained as described in Section II. In contrast to £-Ge, 
but consistent with previous simulations [lcl. . the principal peak in S(k) is strikingly split. 
The calculations are in excellent agreement with experiments carried out on as-quenched a- 
Ge at T = 300 K0|; in particular, the split principal peak seen in experiment is accurately 
reproduced by the simulations. 

We have also calculated a number of other quantities for both £-Ge and a-Ge, including 
pair distribution function g(r), and the electronic density of states n(E). For £-Ge, we 
calculated n(E) using the Monkhorst-Pack mesh with gamma point shifting (one of the 
meshes recommended in the VASP package) . The resulting n(E) is generally similar to that 
found in previous calculations Q, Q] provided that an average is taken over at least 5-10 
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liquid state configurations. Our n(E) for a-Ge [calculated using a shorter averaging time 
than that used below for S(k, u)] is also similar to that found previously^]]. Our calculated 
<?(r)'s for both £-Ge and a-Ge, as given by the VASP program, are similar to those found 
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in Refs. M and Il2j. The calculated number of nearest neighbors in the first shell is 4.18 
for a-Ge measured to the first minimum after the principal peak in g(r). For £-Ge, if we 
count as "nearest neighbors" all those atoms within 3.4Aof the central atom (the larger of 
the cutoffs used in Ref. [h]]) we find approximately 7.2 nearest neighbors, quite close to 
the value of 6.9 obtained in Ref. for that cutoff. Finally, we have recalculated the self- 
diffusion coefficient D(T) for £-Ge at T = 1250 K, from the time derivative of the calculated 

n 

mean-square ionic displacement; we obtain a result very close to that of Ref. [12J . 

B. S(k,u) for £-Ge and a-Ge 

1. i-Ge 

Fig. 3 shows the calculated ratio S(k,u))/S(k) for £-Ge at T = 1250 K, as obtained using 
the averaging procedure described in Section II. We include a resolution function [eqs. (J22J) 
and (J23|) ] of width Tiuj = 2.5 meV, the same as the quoted experimental widthj^lj]. In Fig. 4, 
we show the same ratio, but without the resolution function (i. e., with u = 0). Obviously, 
there is much more statistical noise in this latter case, though the overall features can still 
be distinguished. 

To interpret these results, we first compare the calculated S(k,u) in £-Ge with hydrody- 
namic predictions, which should be appropriate at small k and uo. The This prediction takes 
the form (see, for example, Ref. 5^|): 

S(k,u) 7-1 / 2Drk 2 \ 
n S{k) 7 \u 2 + [D T k 2 ) 2 ) 

1 / Tk 2 Tk 2 \ 

+ 7 V (cu + c s k) 2 + (rp) 2 + (u - c s kf + {Tk 2 } 2 )' ( > 

Here 7 = cp/cy is the ratio of specific heats and constant pressure and constant volume, Dt 
is the thermal diffusivity, c s is the adiabatic sound velocity, and T is the sound attenuation 
constant. and V can in turn be expressed in terms of other quantities. For example, 
Dt = Kr/ipcp), where kt is the thermal conductivity and p is the atomic number density. 
Similarly, r = ~ [0(7 — l)/7 + 6], where a = nr/(pcv) and b is the kinematic longitudinal 
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viscosity (see, for example, Ref. 3|, pp. 264-66). 

Eq. (|24j) indicates that S(k,u) in the hydro dynamic regime should have two propagating 
peaks centered at uj = ±c s k, and a diffusive peak centered at uj = and of width determined 
by Drp- The calculated S(k, uj) / S(k) for the three smallest values of k in Fig. 3, does show the 
propagating peaks. We estimate peak values of huj ~ 10 meV for k = 5.60 nm -1 , huj ~ 11 
meV for k = 7.92 nm -1 , and (somewhat less clearly) huj ~ 13 meV for k = 9.70 nm -1 . The 
value of c s estimated from the lowest k value is c s ~ 2.7 x 10 5 cm/sec. (The largest of these 
three k values may already be outside the hydrodynamic, linear-dispersion regime.) 

These predictions agree reasonably well with the measured S(k, uj) obtained by Hosokawa 
et a/[31], using inelastic X-ray scattering. For example, the measured sound-wave peaks for 
k = ±6 nm -1 occur near 10 meV, while those k = ±12 nm -1 occur at Huj = 17.2 meV, 
Furthermore, the integrated relative strength of our calculated sound-wave peaks, compared 
to that of the central diffusion peak, decreases between k = 7.92 nm -1 and 12.5 nm -1 , 
consistent with both eq. (J2i|) and the change in experimental behavior|31] between k = 6 
nm -1 and 12 nm -1 . 

Because S(k, uj) in Fig. 3 already includes a significant Gaussian smoothing function, a 
quantitatively accurate half-width for the central peak, and hence a reliable predicted value 
for Dt, cannot be extracted. A rough estimate can be made as follows. For the smallest k 
value of 5.6 nm -1 , the full width of the central peak at half-maximum is around 7.5 meV. If 
the only broadening were due to this Gaussian smoothing, the full width would be around 
2hujQ = 5.5 meV. Thus, a rough estimate of the intrinsic full width is ~ V7.5 2 — 5.5 2 = 5 
meV ~ 2hDTk 2 . This estimate seems reasonable from the raw data for S(k,uj) shown in 
Fig. 4. Using this estimate, one obtains Dt ~ 1.3 x 10~ 3 cm 2 /sec. 

The hydrodynamic expression for S(k,uj)/S(k) was originally obtained without consid- 
eration of the electronic degrees of freedom. Since £-Ge is a reasonably good metal, one 
might ask if the various coefficients appearing eq. (J2i|) should be the full coefficients, or 
just the ionic contribution to those coefficients. For example, should the value of Dt which 
determines the central peak width be obtained from the full cp, cy, and kt, or from only the 
ionic contributions to these quantities? For £-Ge, the question is most relevant for kt, since 
the dominant contribution to cp and cy should be the ionic parts, even in a liquid metalQ]. 
However, the principal contribution to kt is expected to be the electronic contribution. 

We have made an order-of-magnitude estimate of Dt using the experimental liquid num- 
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ber density and the value Cp = (5/3)fcn per ion, and obtaining the electronic contribution 
to kt from the Wiedemann- Franz law|5J] together with previously calculated estimates of 



the electronic contribution 



12j . This procedure yields D T ~ 0.1 cm 2 /sec, about two orders 



of magnitude greater than that extracted from Fig. 3, and well outside the possible errors 
in that estimate. We conclude that the D T which should be used in eq. ()24|1 for £-Ge (and 
by inference other liquid metals) is the ionic contribution only. 

In support of this interpretation, we consider what one expects for S(k,u) in a simple 
metal such as Na. In such a metal, ionic motions are quite accurately determined by effective 
pairwise screened ion- ion inter actions [4]. Since the ionic motion is determined by such an 
interaction, the S(k, u) resulting from that motion should not involve the contribution of 
the electron gas to the thermal conductivity. Although £-Ge is not a simple metal, it seems 
plausible that its S(k, us) should be governed by similar effects, at least in the hydrodynamic 
regime. This plausibility argument is supported by our numerical results. 

For k beyond around 12 nm , the hydrodynamic model should start to break down, 
since the dimensionless parameter ujt (where r is the Maxwell viscoelastic relaxation time) 
becomes comparable to unity. At these larger fc's, both our calculated and the measured [3 ll] 
curves of S(k,u)/S(k) continue to exhibit similarities. Most notable is the existence of a 
single, rather narrow peak for k near the principal peak of S(k), followed by a reduction in 
height and broadening of this central peak as k is further increased. This narrowing was first 
predicted by de Gennes|55j]. In our calculations, it shows up in the plot for k = 20.9 nm -1 , 
for which the half width of S(k,u)/S(k) is quite narrow, while at k = 28.5 and 30.7 nm -1 , 
the corresponding plots are somewhat broader and lower. By comparison, the measured 
central peak in S(k, u) / S(k) is narrow at k = 20 nm -1 and especially at k = 24 nm -1 , while 
it is broader and lower at k = 28 nm _1 [31]. 

The likely physics behind the de Gennes narrowing is straightforward. The half-width of 
S(k, oj) is inversely proportional to the lifetime of a density fluctuation of wave number k. 
If that k coincides with the principal peak in the structure factor, a density fluctuation will 
be in phase with the natural wavelength of the liquid structure, and should decay slowly, 
in comparison to density fluctuations at other wavelengths. This is indeed the behavior 
observed both in our simulations and in experiment. 

In further support of this picture, we may attempt to describe these fluctuations by a 
very oversimplified Langevin model. We suppose that the Fourier component p(k, t) [eq. 
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([18)1] is governed by a Langevin equation 

p(k,t) = -£p(k,t)+r)(t). 



(25) 



Here the dot is a time derivative, £ is a constant, and r] (t) is a random time-dependent "force" 
which has ensemble average (rj) = and correlation function {r](i)r]* (t')) = A5(t — t'). Eq. 
fl23|) can be solved by standard methods (see, e. g., Ref. j3] for a related example), with 
the result (for sufficiently large t) 

(p(k,t)p*(k,t + r)) = ^exp(-|r|0. (26) 

According to eq. (fTTj) . S(k, uj) is, to within a constant factor, the frequency Fourier transform 
of this expression, i. e. 

/oo 
(nA/£) exp(iujr) exp(— |t|£)c?t, (27) 
-oo 

or, on carrying out the integral, 



S(k,u) oc-^--. (28) 



7rA 

This is a Gaussian function centered at uj = and of half-width £. On the other hand, the 
static structure factor 

71 A 

S{k) oc Lim T ^ (p(k, t)p*{k, t + r)) = — . (29) 

Thus, if the constant A is independent of k, the half-width £ of the function S(k, uj) at wave 
number k is inversely proportional to the static structure S(k). This prediction is consistent 
with the "de Gennes narrowing" seen in our simulations and in experiment [sil] . 

To summarize, there is overall a striking similarity in the shapes of the experimental and 
calculated curves for S(k,uj)/ S(k) both in the hydrodynamic regime and at larger values of 
k. 



2. a-Ge 



We have also calculated the dynamic structure for our sample of a-Ge at T = 300K. 
The results for the ratio S(k,uj)/S(k) are shown in Figs. 5 and 6 for a range of k values, 
and, over a broader range of uj, in Fig. 7. Once again, both S(k, uj) and S(k) are averaged 
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over different values of k of the same length, as described above. We have incorporated a 
resolution function of width Tiojq = 2 meV into S(k,u). This width is a rough estimate for 
the experimental resolution function in the measurements of Maley et a/|56j; we assume 
it to be smaller than the liquid case because the measured width of the central peak in 
S(k, u) I S(k) for a-Ge is quite small. 

Ideally, our calculated S(k, u>) should be compared to the measured one. However, the 
published measured quantity is not S(k,u) but is, instead, based on a modified dynamical 
structure factor, denoted G(k,u), and related to S(k,u) byjsfj 

G(k, lj) = (%) r) S(k, u>). (30) 



,k 2 J \n(w,T) + 1 

Here C is a k- and (^-independent constant, and n(cu,T) = \j\^ l ^l k B T _ 1] j s the phonon 
occupation number for phonons of energy E = Tiu at temperature T. The quantity plotted 
by Maley et a/[^f| is an average of G(k,u) over a range of k values from 40 to 70 nm" 1 . 
These workers assume that this average is proportional to the vibrational density of states 
n V ib(u). The measured n v ib(uj) as obtained in this way 56] is shown in Fig. 8 for two different 
amorphous structures, corresponding to two different methods of preparation and having 
differing degrees of disorder. 

In order to compare our calculated S(k, uj) to experiment, we use eq. (|3T?j) to infer G(k, u), 
then average over a suitable range of k. However, in using eq. ()30|) . we use the classical form 
of the occupation factor, n{Tiuj) + 1 ~ A^T /hu, This choice is justified because we have 
calculated S(k, u) using classical equations of motion for the ions. We thus obtain for the 
calculated vibrational density of states 

n vib (uj) w (j^) fel S(k,u). (31) 

In Fig. 8 we show two such calculated plots of n vib (u), as obtained by averaging eq. (}3Tj) 
over two separate groups of k's, as indicated in the caption 5^. For comparison, we also 
show n„ifc(co>) for a-Ge as calculated in Ref. directly from the Fourier transform of the 
velocity-velocity autocorrelation function. 

The calculated plots for n v n,{oj) in Fig. 8 have some distinct structure, which arises from 
some corresponding high frequency structure in S(k,u). The plot of n V ib(uj) for the group 
of smaller k's has two distinct peaks, near 8 meV and 29 meV, separated by a broad dip 
with a minimum near 18 meV. The plot corresponding to the group of larger k's has similar 
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structure and width, but the dip is less pronounced. The two experimental plots also have 
two peaks separated by a clear dip. The two maxima are found around 10 and 35 meV, while 
the principal dip occurs near 16 meV. In addition, the overall width of the two densities of 
states is quite similar. 

The reasonable agreement between the calculated and measured n vib (u) suggests that our 
ab initio calculation of S(k, u) for a-Ge is reasonably accurate. The noticeable differences 
probably arise from several factors. First, there are several approximations involved in going 
from the calculated and measured S(k,u>ys to the corresponding ri^(a;)'s, and these may 
be responsible for some of the discrepancies. Secondly, there may actually be differences 
between the particular amorphous structures studied in the experiments, and the quenched, 
then relaxed structure considered in the present calculations. (However, the similarities 
in the static structure factors suggest that these differences are not vast.) Finally, our 
calculations are carried out over relatively short times, using relatively few atoms; thus, 
finite-size and finite-time effects are likely to produce some additional errors. Considering 
all these factors, agreement between calculation and experiment is quite reasonable. 

Previous ab initio calculations for a-Ce Ujf have also obtained a vibrational density of 
states, but this is computed directly from the ionic velocity-velocity autocorrelation function 
rather than from the procedure described here. The calculations in Ref. [lOj do not require 
computing S(k,u). In the present work, by contrast, we start from our calculated S(k,u), 
and we work backwards to get n v ib(uj). In principle, our S(k,u) includes all anharmonic 
effects on the vibrational spectrum of a-Ge, though in extracting n V ib(u) we assume that 
the lattice vibrates harmonically about the metastable atomic positions. In Fig. 8, we also 
show the results of Ref. for n vi b(uj) as obtained from this correlation function. They are 
quite similar to those obtained in the present work, but have a somewhat deeper minimum 
between the two principal peaks. 

The quantity n V ib(uj) could, of course, also be calculated directly from the force constant 
matrix, obtained by assuming that the quenched configuration is a local energy minimum 
and calculating the potential energy for small positional deviations from that minimum using 
ab initio molecular dynamics. This procedure has been followed for a-GeSe2, for example, by 
Cappelletti et al[58]. These workers have then obtained S(q,u) versus q from their n V ib{oj) 
at selected values of u, within a one-phonon approximation. However, as noted above, the 
present work produces the full S(k, u) and thus has, in principle, more information than 
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n vib {u). 



V. DISCUSSION AND CONCLUSIONS 



The results reviewed here show that ab initio molecular dynamics can be used to calcu- 
late the dynamic structure factor S(k,u) for both liquid and amorphous semiconductors. 
Although the accuracy of the calculated S(k, uj) is lower than that attained for static quan- 
tities, such as S(k), nonetheless it is sufficient for comparison to most experimental features. 
This is true even though our calculations are limited to 64-atom samples and fewer than 20 
ps of elapsed real time. 

We have presented evidence that the calculated S(k,u)/S(k) in £-Ge agrees qualitatively 



with measured by inelastic X-ray scattering 3l|, and that the one calculated for a-Ge leads 
to a vibrational density of states qualitatively similar to the quoted experimental one 
Since such calculations are thus shown to be feasible, the work reviewed here should spur 
further numerical studies, with longer runs on larger samples, to obtain even more detailed 
information. Furthermore, we can use these dynamical simulations to probe the underly- 
ing processes at the atomic scale which give rise to specific features in the measured and 
calculated S(k,u). 
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FIG. 1: Static structure factor S(k) for l-Ge at T = 1250K, just above the experimental melting 
temperature. Full curve: present work, as calculated using the generalized gradient approximation 
(GGA; see text). Dashed curve: present work, but using the local density approximation (LDA; 
see text). Open circles: measured S(k) near T = 1250 K, as given in Ref. 51]. 
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FIG. 2: Full curve: Calculated S(k) for a-Ge at T = 300K, as obtained using the GGA. Structure 
is prepared as described in the text. Open circles: measured S(k) for a-Ge at T = 300-fT, as given 
in Ref. 
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FIG. 3: Calculated ratio of dynamic structure factor S(k, uo) to static structure factor S(k) for £-Ge 
at T = 1250-ftT for several values of k, plotted as a function of u, calculated using ab initio molecular 
dynamics with a MD time step of 10 fs. For clarity, each curve has been vertically displaced by 
0.05 units from the curve below. In each case, the plotted curve is obtained by averaging both 
the calculated S(k, u) and the calculated S(k) over all values of k of the same length. We also 
incorporate a Gaussian resolution function of half-width Hujq = 2.5 meV, as in eqs. ()22l) and (|23|) . 
This value of loq is chosen to equal the quoted experimental resolution |3J| . 
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FIG. 4: Same as in Fig. 3, but without the resolution function. 
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FIG. 5: Calculated S(k,u)/S(k) for a-Ge at T = 300 K at k < 35 nn -1 , plotted function of 
u) Again, each curve has been vertically displaced by appropriate amounts from the one below it, 
as evident from the Figure, and both S(k, to) and S(k) have been plotted after an average over all 
k's of the same length. We also incorporate a Gaussian resolution function of half-width Tiloq = 2 



meV. This value is chosen to give the best results for n v ib(ui) as measured by Ref. 
step here is 10 fs. 
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FIG. 7: Calculated S(k, u)/S(k) as in Figs. 5 and 6 but including higher frequencies u. At huj = 60 
meV, the curves are arranged vertically in order of increasing frequency. Each curve is vertically 
displaced by 0.0005 units from the one below it. 
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FIG. 8: Full curve: calculated vibrational density of states n v n,{ui), in units of 10 -3 states/meV. 
nvib(u) is obtained from the resolution-broadened S(k,u) and S(k) of the previous Figure, using 
the formula n v n,{io) = (G ca i c (k,u>)) , where G ca i c (k,uj) is given by eq. (16), and the averaging is 
carried out over the three magnitudes of k near 40 nm" 1 for which we have computed S(k,uj). 
Dashed curve: same as full curve, but calculated by averaging over the six magnitudes of k near 90 
nm" 1 for which we have computed S(k,co). The open circles and open stars denote the measured 
n v ib(uj), as reported in Ref. |^6j for two forms of a-Ge. Finally, the open diamonds denote n v ib(u>) as 
calculated in Ref. [l£J] from the ionic velocity- velocity autocorrelation function (dot-dashed curves). 
In all plots except that of Ref. [ljj, n v ib(io) is normalized so that J^ max n{uj)d{%uj) = 1. uj max is 
the frequency at which n v ib(u) — > 0, and is estimated from this Figure by extrapolating the right 
hand parts of the solid and dashed curves linearly to zero. 
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